function [elements,processInfo] = InitializeSolutionStep...
    (elements,model,processInfo)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    xi = e.getGaussPoint(0);
    w = e.getGaussWeight(0);
    jac = f.jacXAt([xi(1) xi(2)]);
    detJac = jac(1)*jac(4) - jac(2)*jac(3);
    B = zeros(3,6);
    B = e.computeB(f,B,[xi(1) xi(2)]);
    C = elements(i).constitutiveMatrix;
%     %compute Xi
%     if elements(i).isCrack == 1
%         due = ue - elements(i).ue_old;
%         dxi = elements(i).Kxixi\(elements(i).renh-...
%             elements(i).Kxid*due);
%         elements(i).Xi = elements(i).Xi_old+dxi;
%     end
    %compute sigma
    if elements(i).isCrack==0
        elements(i).stress = C*B*ue;
    else
        epsilon = B*ue + elements(i).Gc*elements(i).J*elements(i).Xi;
        elements(i).stress = C*epsilon;
    end
    %update ue
    elements(i).ue_old = ue;
    %%
    %compute renh at element's crack
    if elements(i).isCrack==1
        %compute TGamma and update Kxixi
        lGamma = elements(i).lGamma;
        J = elements(i).J;
%         elements(i).Xi
        u_mu = J*elements(i).Xi;
        t = elements(i).ft+elements(i).H*u_mu(1);
        TGamma = [max(t,0);0];
%         TGamma = max(t,0)*elements(i).crack_normal;
        if t < 0
            elements(i).Kxixi = w*elements(i).Ge'*...
                C*elements(i).Gc*abs(detJac);
        else
            Ctilde = [elements(i).H 0;0 0];
            elements(i).Kxixi = w*elements(i).Ge'*...
                C*elements(i).Gc*abs(detJac)+...
                lGamma*J'*Ctilde*J;
        end
        
        %compute renh
        renh = -lGamma*elements(i).J*TGamma;
        renh = renh - w*elements(i).Ge'*...
            elements(i).stress*abs(detJac);
        elements(i).renh = renh;
%         renh
    end
end
%elements = calculateOnsetCriteria(elements,model);
end
